Majorana fermions in one-dimensional spin-orbit coupled Fermi gases 
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We theoretically study trapped one-dimensional Fermi gases in the presence of spin-orbit coupling 
induced by Raman lasers. The gas changes from a conventional (non-topological) superfluid to a 
topological superfluid as one increases the intensity of the Raman lasers above a critical chemical- 
potential dependent value. Solving the Bogoliubov-de Gennes equations self-consistently, we calcu- 
late the density of states in real and momentum space at finite temperatures. We study Majorana 
fermions (MFs) which appear at the boundaries between topologically trivial and topologically non- 
trivial regions. We linearize the trap near the location of a MF, finding an analytic expression for 
the localized MF wavefunction and the gap between the MF state and other edge states. 
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I. INTRODUCTION 

Majorana fermions (MFs) , exotic excitations which are 
their own antiparticles, have attracted a great deal of at- 
tention recently [1]. Condensed matter systems with MFs 
possess degeneracies that are intrinsically nonlocal, and 
can be manipulated geometrically. They can, in prin- 
ciple, be used to make a robust quantum computer [2]. 
Condensed matter theorists have proposed various ways 
to explore MFs during the past several years [3-8] . Four 
experimental groups have recently reported evidence of 
MFs in semiconducting wires on superconducting sub- 
strates [9], In those experiments, spin-orbit (SO) cou- 
pling was important. Here we study MFs in a related 
cold atom system. 

Two groups [11, 12] have successfully generated SO 
coupled Fermi gases based on a Raman technique pio- 
neered by Spielman's group at NIST [10]. Several the- 
oretical groups have proposals for creating and probing 
MFs in these SO coupled Fermi gases [13-16]. We build 
upon the studies of Jiang et al. [14] and Liu et al. [15], 
which find MFs in a ID geometry. 

We study a ID (pseudo) spin-1/2 Fermi gas with point 
interactions. In the presence of Raman lasers, the energy 
spectrum has two helical bands. We study this two-band 
model in a harmonic trap. Solving the Bogoliubov-dc 
Gennes (BdG) equations self-consistently, we calculate 
the density of states (DOS) in real and momentum space 
at finite temperatures. We linearize the trap near the 
location of a MF, finding an analytic expression for the 
localized MF wavefunction and the gap between the MF 
state and other edge states. 

Our numerical calculations extend the similar studies 
of Ref. [15]. We explore a larger range of temperatures, 
and delve deeper into the physics near the MFs. We also 
investigate a truncated one-band model. 

One concern with mean-field calculations such as ours, 
is that they are unable to capture the large phase fluc- 
tuations found in ID. As shown by Ref. [19], the MF 
physics is robust against these fluctuations. Moreover, 



an actual experiment would be performed on a bundle of 
weakly coupled tubes [17]. This latter setting also avoids 
issues of number conservation [19]. Our ID model faith- 
fully describes the properties of a single tube within such 
a bundle when the tunneling is weak. 

This paper is organized as follows. In Sec. II, we dis- 
cuss the homogeneous gas: We start with the two-band 
model, and in Sec. 11(A) show how it relates to a one- 
band model with p-wave interactions. In Sec. 11(B), we 
describe the band structure and topology of the two-band 
model. In Sec. Ill, we calculate the properties of trapped 
gases: In Sec. III(A), we write the BdG equations and 
self-consistently calculate the order parameter and den- 
sity. In Sec. III(B), we visualize the MFs by calculating 
the DOS in real space and momentum space. In Sec. 
III(C), we introduce MF operators and construct the lo- 
calized MF states. In Sec. III(D), we linearize the trap 
near the location of a MF, finding an analytic expression 
for the localized MF wavefunction and the gap between 
the MF state and other edge states. Finally we conclude 
in Sec. IV. 



II. HOMOGENEOUS GAS 

We start from the Hamiltonian of the ID (pseudo) 
spin-1/2 Fermi gases with chemical potential fx, 



H = 



W(x)(H Q (x) -/i)tf(x) )dx + H 



(1) 



where ^(x) = (^)-j-(x), xj}±(x)) T annihilates the spin-up 
and spin-down states. In an experiment, ip^ and ip± cor- 
respond to two different hyperfine states of a fermionic 
atom such as 4Q K. The single-particle Hamiltonian 



H (x) 



-d 2 

2m x 



d x a. 



Ml. 



E r can be en- 



gineered by Raman lasers [10] , whose intensity is charac- 
terized by the Rabi frequency f2. The recoil momentum 



of the Raman lasers is hkt, E r 



2m 



is the recoil en- 



ergy, and cr = (a x , cr y ,cr z ) is the vector of Pauli matrices. 
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For ultra-cold fermions, the interaction may be modeled 
by Hj = guj J tpttpyip^ip^dx, with coupling constant g^D- 
This coefficient can be related to the three-dimensional 
scattering length and the geometry of the confinement 
[20]. In a typical experiment, \gm\ ~ 70aoE r [17], where 
ao is the Bohr radius. We restrict ourselves to attractive 
interactions, g±o < 0. We note that if we rotate our spin 
basis (a x — > u z ,cr z — > a y ) and identify Z = ^ as a Zee- 

man field, and a = 2 ^ as the SO coupling strength, we 
recover the Hamiltonian of a semiconducting wire. Note 
Hi is very different for a wire [5]. In the following sec- 
tions we explore the physics of Eq. (1). 



A. One-band model 

To get insight into Eq. (1), we first consider an ap- 
proximation where we truncate to a single band. We 
emphasize however that in all other sections, we work 
with the full two-band Hamiltonian. 

The physics of the single particle Hamiltonian 
is most transparent in momentum space, H = 

E fc *I i^+Er- + *fc, where * fe = 

(V'fefj "0h) T - This Hamiltonian is readily diagonalized by 



Ck 



cos^ — sin^ 
sin^ cos^ 



( ^fct 



(2) 



with 



E fc (W^ + E+d\d h 



H 



2WkZ> yielding 
ik j . The energy spectrum 
has two helical bands, illustrated in Fig. 1. If the 
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FIG. 1: (color online) Band structure of a ID (pseudo) spin- 
1/2 gas. The red (dashed) curves are the bare bands in the 
absence of SO coupling. The blue (thick) curves are the upper 
band E+ and lower band E- in the presence of SO coupling, 
with the coupling strength M~l/E r — 1. 

effective chemical potential p, = /i — E r <C 4P , only the 
lower band E_ is filled with fermions. Projecting the 
interactions into this band, we find 



H} B = giD ( f *<j 



,t 



kqq' 



V k q>Ck q ,Ck 



(3) 



where g±n = 9id/Lid, with Lid the length of the gas. 
The fermionic anti-commutation relation, c 



...t 



„t 



■+q 



implies that the interaction coefficient Vkq 



f-9 f+9 

is odd with respective to q, Vk q 
At zero center of mass momentum, V q 



1 • "l/^jV^t/z- 
2 b 1 2 



Vk=0,q 



— \ In Fig. 2, we plot V4 as a function 

of g. The dependence on k is weak for k < k^. 

The interaction in Eq. (3) is separable. Given that 
gio < 0, this interaction can lead to pairing with 
zero center of mass and an order parameter A q = 

Ti( e - H/k i> T ...) . ^ , , 

-q'/, wiicic _ Tr ^ e _j?/fc ( ,Ts is tne 
thermal average, fc;, is the Boltzman constant and T is 
the temperature. The mean-field interaction becomes 

By virtue of the symmetry of V q , the order parameter 



giDV q Y^ q >(Vq'C-q'Cq')i wherC 



E 9 ( A <j4 C -9 + A q C -1 C l) ~91d\ Eg^g( C -« 



has a p-wave symmetry A_ 



-A„. As is well estab- 



lished, such a p-wave superfluid may possess Majorana 
edge modes [4]. We will discuss these Majorana modes 
at length in the two-band model. 
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FIG. 2: (color online) Interaction coefficient Vkq versus di- 
mensionless momentum g/fcz, for hQ/E r = 2. The blue 
(thick), green (dashed) and red (dotted) curves correspond 
to k = 0,0. 5&l and fez, respectively. 



B. Two-band model 

While the one-band model connects the SO cou- 
pled gases and p-wave superconductors, we will focus 
on the richer two-band model in the remainder of the 
manuscript. Within the mean-field approach, the inter- 
action term is bilinear 



H, 



9lD 



A 



(x) U>\i/>1 + 



A(x)< 
giD 



(4) 



dx, (5) 



where the order parameter A(x) = g\r> (VM-^t) i g as ~ 
sumed to be real. Defining the operator W (x) = 



3 



ytp^(x), ip^(x), tpi(x), ipf(x)J , the Hamiltonian can be 
written as, 



±E ± /E r 



±E t /E r 



H = 



where 



^¥(x)HV{x) 



A(xf 

9lD 



dx + \ (T_ 



+ T+) ,(6) 



u 



T+ 



2m 



HQ 
Tr 



2m x 



- M 
A(x)t x <j z , 



ih 2 k L 



'-d x r z a z 



(8) 



The Pauli matrices er, r operate in the spin subspace and 
particle-hole subspace respectively, 

















... 











/I 

V 



-1/ 



-1 



(9) 



(10) 



-I J 



The elementary excitations can be found by solving the 
BdG equations UW = EW . When A(x) = A is spatially 
homogeneous, one can write the BdG equations in mo- 
mentum space as HkW(k) = E(k)W(k), where Ht is the 
4x4 matrix produced by replacing — id x — > k in Eq. (7). 
The excitation spectrum E(k) is most simply calculated 
by squaring Tik twice, and extracting the characteristic 
polynomial [7]. This procedure yields 



= £q + 2E r h k /m + H CI /4 + A 



± hJ8E r e^k 2 /m + fl 2 ei + fl 2 A 



(11) 



where en = 4r^- — U ■ The four bands 

E + (k),E_{k),-E_(k),-E + {k), as shown in Fig. 
3, correspond to the four eigenvectors W p+ (k), W p ~(k), 
W h -(k) 7 W h+ {k). The Hamiltonian U has the intrinsic 
symmetry, {%,T y } = 1. Given two eigenvectors W p± 
with eigenvalues E ± , one can always construct the 
other two W h± = ir y W p± with eigenvalues —E ± . We 
therefore denote, 



W p+ (k) = 


( u kf 


u tl 




4 t y 


(12) 


w p -(k) = 


( u feT 


u ki 






(13) 


W h ~(k) = 


( v kV 






) T 


(14) 


W h+ (k) = 


( v k V 




+ + 

U fcf U kl 


) T - 


(15) 



The unitary condition on the 4 by 4 matrix 
(W p+ (k), W p -{k),W h -{k) 7 W h+ (k)) also leads to the 



equalities = (ut^)* and v± = -(vt^)* 



\ (a) V 
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FIG. 3: (color online) Band structure of homogeneous 
gas. From the top to the bottom, the four bands are 
E+,E-,—E-,—E+ respectively. The parameters are fi = 
E r , A = 0.5£ r , and (a)fi = 0, (b)fi = 0.5£ r , (c)fi = E r , 
(d)fi = 1.5E r . 



In Fig. 3, the spectrum is shown for a range of pa- 
rameters. One important feature is the k = gap Eq = 
2E_(k = 0) = 2\G\, where G = ^jj? + A 2 -Ml/2. When 
f2 = 0, the two positive energy bands touch at k = 0, and 
Eq > 0. The gas is in the same universality class as a con- 
ventional s-wave superconductor. Increasing the Raman 
laser strength such that < Ml < 2^J A 2 + j} 2 separates 

the two bands and reduce Eq. At Ml = 2y/ A' 2 + /i 2 , 
Eq is zero, and there is a topological transition. Once 
Ml > 2y/ A 2 + p, 2 , Eq is again positive, but the gas is 
no longer a conventional supcrfluid, instead it has a non- 
trivial topological invariant. 

The relevant topological invariant is a Berry phase. 
Eqs. (12-15) can be thought of as maps from the real line 
(— oo < k < oo) to the space of unit vectors in SU(4). 
One can generate a closed path by taking 

C: W p +{-oo) -> W p+ {oo) = W p -(-oo) -> 

W p "(+oo) =W p+ {-oo). (16) 

The equalities follow from noting that up to phases 

W p+ (-oo) = W p -(oo) = (0,1,0, 0) T (17) 
W p -(~oo) = W p+ (oo) = (1,0,0, 0) T . (18) 

Given this closed path, one can dehne the Berry phase 

(19) 

= i I (W p+ (k))* -d k W p+ (k)dk 

J — oo 

/oo 
(W p -{k))* ■ d k W p -(k)dk. (20) 
- oo 

In the case of a gauge which is not smooth, one would 



i J> W* ■ dkWdk 
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instead use 

(oo 
(W p+ (k))* -W p+ {k~8k) 
k— — oo 

OO \ 

x (W p -{k))* ■ W p -(k-8k)\ (21) 



k=— oo 



x (W p+ (-oo))* ■ W p -(oo)(W p -(-oo))* ■ W p+ (oo). 

Since H has real valued matrix elements, e %1 = ±1. The 
Berry phase 7 is only well defined if the spectrum has 
no degeneracies. We restrict < 7 < 2ir. For G > 0, 
we find 7 = 7r. For G < 0, 7 = 0. Somewhat counter- 
intuitively the 7 = sector corresponds to the "topologi- 
cally non-trivial" state analogous to a ID spinless p-wave 
superconductor . 



III. TRAPS 

In this section we will solve the BdG equations for a 
trapped gas. The qualitative features of our results can 
be anticipated by treating the system as locally homo- 
geneous: the properties at position x will be reminis- 
cent of those corresponding to a homogeneous gas with 
chemical potential fl(x) = fb — V(x). Within this local 
density approximation (LDA), one can define a function 
G{x) = yjp,{x) + A(x) 2 - Ml/2, where G(x) = corre- 
sponds to the boundaries between topologically distinct 
regions. One expects there will be Majorana excitations 
at the boundaries. We will use numerical solution of the 
BdG equation to explore this physics beyond the LDA. 
Further, in Sec. HI(D) we will linearize the BdG equa- 
tions about the points G{x) = 0, and analytically inves- 
tigate these Majorana modes, without making a LDA. 



A. Order parameter and density 

In the presence of a trap, we write the BdG equations 
in real space, 



where 



UtrapW n {x) = E„W„{X), 



2m 



(22) 



dt - A + V{X) T; 



ih 2 kL „ Ml 

+ o x T z a z + —T z a x + A(x)T x a z . (23) 

m 2 

The eigenvectors W n {x) come in pairs, W p {x) and 
W£(x), which correspond to eigenvalues E n > and 

— E n , 

W ni x ) = (Unt(x),U nl (x),V nl (x),V n1 -(x)) T (24) 

W*{x) = (v; it (x),v* ni (x),u* ni (x),u* nt (x)) T . (25) 



To make contact with our previous discussion, we note 
that in the spatially homogeneous case, n can be repre- 
sented by a momentum k n and a sign e„ = ±, so that 
W p (x) = 



%X {<V<lv V knV V l^y and W n( X ) 

e^((^ t )*,«- 4 r,(u^)*,(«^ t )*) T . One then re- 
covers Eqs. (12)-(15). 

Fixing {/2, fl, <7i£>}, we solve Eqs. (22) itcrativcly. We 
discretize space into n gr id equally spaced points, and 
use a finite difference method with a pseudo-spectral 
scheme to represent Ht ra p as a 4?i gr id by 4n gr ;d ma- 
trix. In the jth iteration, we numerically diagonalizc 
the matrix W^} ap with the order parameter A^\x). We 
start with a constant A^(x) = A . We extract the 
eigenvectors Wn\x) and calculate the order parameter 

AW+D(x) = 9inE p (u^v:f{^a) + v n fu^n)), 
where is the annihilation operator of the Bogoliubov 
particle. At temperature T, = l/(e B «/ kt ' T + 1). 

Then we diagonalize Ti^ra^p an d repeat the procedure. 
We stop iterating when \A^ +1 ^(x) — A^(x)\ falls be- 
low a threshold. The final convergent order parameter 
i'^Ji) is largely independent of n gr id and A(°)(a;) when 
"grid > 1200. In the Appendix we explore the conver- 
gence with the real space grid size ro gr id- 
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FIG. 4: (color online) Profiles of order parameter A(a;) = 
Sid J2„{ u ni(x)v^(x){^i) + v'^UnYWdUn)) (dashed 
curves) and density n(x) = J2 na (I*w(a0| 3 (£n£i) + 
\una {x)\ 2 {£}„£,n)) (solid curves) at temperatures T = 
0, 0.1i? r , 0.2i? r , 0.3-EV. Other parameters are gm = 
-0.03E r L, htt = 2E r ,\ = 4, k L L = 100, (a)/i = E r ,(b)fx = 
2.5E r . 

The order parameters and density profiles for a gas 
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in a harmonic trap V(x) = X(x/L) 2 E r are shown in 
Fig. 4, where the dimensionless parameter A = 4 char- 
acterizes the stiffness of the trap, and 2L is the simula- 
tion length with fc^L = 100. We choose the Rabi fre- 
quency to be Ml = 2E r , and take gio = —0.03E r L, 
corresponding to the dimensionless interaction strength 
P = m\giD\/ti 2 nQ ~ 2, where no is the central density at 
zero temperature. For comparisons, experiments on ID 
Fermi gases at Rice have j3 ~ 3 [17]. If E r /H — 50kHz 
(a typical experimental value), then these parameters 
correspond to a trap with small oscillation frequency 
lu = 2kHz. The order parameter has qualitatively dif- 
ferent behavior if the center of the trap has one or two 
bands occupied. For relatively small chemical potential 
E r - y/(m/2) 2 - A(x) 2 < n < E r + y/(m/2) 2 - A(x) 2 , 
the center of the trap will be topologically non-trivial 
while the wings will be trivial. This regime is illus- 
trated in Fig. 4(a). The order parameter grows near 
the edge of the cloud. This is a feature of ID where, 
due to the divergence of the low energy density of state, 
the interactions are stronger for lower density [21]. For 
p, > E r + a/ (Ml/2) 2 — A(x) 2 , the center will be topolog- 
ically trivial, but there will be a band of the non-trivial 
state further out. Here the order parameter profile is 
quite rich, with a central plateau, surrounded by two 
valleys and two peaks. The central plateau roughly cor- 
responds to where both bands are occupied. The order 
parameter is sensitive to temperature. The bulk A is sig- 
nificantly suppressed and vanishes for T > Q.2E r , The 
density has no notable structure and is nearly indepen- 
dent of temperature for T < 0.3E r . 

B. Density of states (DOS) 

The elementary excitations are encoded in the 
single particle Green function G aa > (x, t, x' , t') = 
i(Ttp a (x, t)i/j a r(x' ,t')) and the associated spectral den- 
sity A a(T ,(x,x',E) = 21m / e lEt G aiJ ,(x,t,x',0), where f 
is the time-ordering operator. A local tunneling experi- 
ment can measure the density of states (DOS) p(E, x) = 
A^(x,x,E) + An(x,x,E). This quantity gives the num- 
ber of single particle states with energy E at position x. 
It can be understood as an application of Fermi's Golden 
rule to the response to a tunneling probe. Within our 
mean-field theory 

p(E,x) = J2 (Pa(-E,x)+(? a (E,x)), (26) 

<7=f4 

where 

p h a (E,x) = J2\vn*(x)\ 2 5(E n -E) (27) 

n 

(%(E,x) = Y,\ u nAx)\ 2 6(E n -E). (28) 

n 

We can similarly introduce a momentum resolved DOS 
p{E,k) = J e ik ^ x - x '\A n {x,x\E) + A u (x,x',E))dxdx', 
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FIG. 5: (color online) Density of states (DOS) in real space 
(left panel) and momentum space (right panel) at T = 
0, 0.lEr, 0.2E r , 0.3E r from the top to the bottom, with or- 
der parameters identical to those in Fig. 4(a). The grey 
(dashed) curves in the left panels is plotted with G(x) — 
\/p(x) 2 + A(x) 2 - m/2, where the zero points of G(x) pin- 
point the position of MFs. The brighter color corresponds to 
the higher spectral weight. 



which can be measured with momentum resolved radio- 
frequency spectroscopy [12]. In the present case 



p(E,k) = {pi(-E,k)+pP(E,k)), (29) 

<r=t4 
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where 



P h AE,k) 



E 

n 

E 



v na {x)e x dx 



Unaix)^ 1 dx 



6{E n -E) (30) 



8{E n -E). (31) 



Finally for comparison, we plot the DOS within a LDA. 
As illustrated in Fig. 7, the LDA prediction for the DOS 
is qualitatively similar to the BdG result. The main dif- 
ference is that the LDA misses physics related to quan- 
tization. In particular, the zero energy modes are not 
spectrally isolated in the LDA. They are, however, still 
located at roughly the same place in space. 






FIG. 6: (color online) Density of states (DOS) in real space 
(left panel) and momentum space (right panel) at T = and 
0.3E r , with parameters identical to those in Fig. 4(b). The 
brighter color corresponds to the higher spectral weight. 

In Figs. 5 and 6 we plot the DOS for the trapped 
gas with the order parameters calculated in Sec. Ill A. 
We also show a dashed curve corresponding to G(x) = 
\Jjl{x) 2 + A(a;) 2 - Ml/2. The point where G(x) = rep- 
resents the boundary between topologically distinct re- 
gions defined in Sec. II B. For the parameters in Fig. 5, 
G(x) — at two locations, and we find that the BdG 
equations have two zero-energy modes, localized near 
these points. As will be discussed later, these modes 
may be interpreted as MFs. They are clearly spectrally 
separated from all other states. Fig. 6 shows the case 
where G(x) = at four locations, representing four MFs. 
The right panels of Figs. 5 and 6 show the momentum 
space DOS. The MF modes sit in a large gap at k = 0. 

As we have shown in Sec. Ill A, the order parameter 
decreases with temperature. In real space, the bulk A 
becomes very small at T = 0.2E r , while A at the edges 
remains large: the MFs at the edges are very clear for 
T < 0.2E r . At T = Q.3E r , the order parameter is nearly 
zero and the gas becomes normal. 

The evolution of the momentum space DOS parallels 
the real space DOS. As temperature is increased from 
T = 0, the gaps at large k shrink. The gap at k = 
remains robust until T = 0.3E r . 



FIG. 7: (color online) Density of states (DOS) at zero tem- 
perature under the local density approximation (LDA). The 
parameters are identical to those in Fig. 4, except A (a;) is 
calculated within the LDA. The brighter color corresponds to 
the higher spectral weight. 



C. Majorana fermions (MFs) 

Here we explore the structure of the zero-energy states 
seen in Fig. 5. From our numerical solutions to the BdG 
equations, we have two wavefunctions 

Wg(x) = e l ^(uo t (x),u ol (x),v ol (x),v ot (x)y (32) 
W h (x) = e^(^ t (x),^( a ;),u* i (.T),^ t ( a :)) T ,(33) 

which define operators 

Co = j (34) 

eJ = eM-^ J («(x)) f • §(*)) dx, (35) 

and obey HtrapWg (x) « UtrapW^x) w 0. The phases 
ipi and ip' 2 are not unique, and the factor in Eq. (32) 
must be introduced to make £g conjugate to £o- By con- 
struction these are fermionic operators {£oj£q} = 1- 
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As zero-energy solutions to the BdG equations, both 
£o and £j commute with H. Hence the ground state 
is degenerate: £ |GSi) = and \GS 2 ) = £ \GS 1 ). These 
two degenerate states can be used as a qubit for quantum 
information processing [2]. 

The operator £g which couples \GS\) to \GS 2 ) is intrin- 
sically nonlocal, with weight at two spatially separated 
points. One can however introduce operators 



Xo 



V2 



Xo = ± ^ lV (& - e- 2iv e5) = J fl(?) ■ *(x)dx(37) 
where 

f (x) = l= e ^(W^x) + e-^W h (x)) (38) 

fo(x) = ±-i-e^«(x)-e-^<(x)), (39) 

with arbitrary phases tp\ = Lp — ip^ and tp 2 = 2<y9 — ^'i + V^- 
By choosing the appropriate (p 2 , these operators can be 
made local. In particular if G(x) = at x = x\, x 2l then 
fo{x) can be chosen to be nonzero only near x\ , and fo{x) 
only near x 2 . 

The operators \o and Xo obey the Majorana algebra: 

Xo = Xo, xt = Xo, {Xo,Xo} = 0, {xo,Xo} = {xo,Xo} = 1- 
They commute with the Hamiltonian. 

Note, as we will use in the next subsection, fo{x) = 
e lip f (iif^(x),Ufi(x),Vfi(x),Vf^(x)) obeys the BdG equa- 
tions, but the resulting Bogoliubov transformation is not 
unitary as it changes the commutation relations. Since 
Xo = Xo> we nave u fcr( x ) = e~~ 2 ' l(p fv*f a (x). For smaller 
systems, coupling between these modes push them away 
from Eq = 0. 



D. Eigen-energies of excited states near a MF 

As seen in Figs. 5-6, the MFs are localized in real 
space and momentum space. Thus we can calculate their 
properties by linearizing the trap around their locations 
in position space, and linearizing momentum around k = 
0. As previously discussed, the locations of the MFs can 
be found via the LDA. There are generally four MFs, 
localized at x m — ±Ly / R± m / XE r , where R± m = ft ± 
^ft 2 fi 2 /4- A r 2 „, with A m = A (a; = x m ). We restrict 
ourselves to the location of one MF, x m = L\J R+ m / XE r . 
We write the linearized BdG Hamiltonian as the sum of 
two terms %u n = Ho + Hi, 

nn 



Uo = — r z a x + A m r x + y/mPjI^^T x cr z (40) 



Hi = X{x - x m )r z - KT z a z , 



(41) 



where A = 2Xx m E r / L 2 and k = T^k^klm. The "inter- 
action" term Hi can be treated as a perturbation, and 



it vanishes as x — ¥ x mi k — s- 0. In the absence of per- 
turbations, %ii n = Ho has two degenerate zero-energy 
states 



y/2, . 

D\ = -— (sin</>, — cos</>, — cos</>, sin<; 



V2, 
2 



T>2 = — (— cos0, sin</>, — sin</>, cost; 



(42) 
(43) 



where sin0 = y/ {Ml/2 + A m )/ftf2. Following the stan- 
dard approach to first-order degenerate perturbation the- 
ory, we diagonalizc the Hamiltonian projected into the 
subspace {T>i,T> 2 }, 



H 



Pi 



H H JV X ,V 2 ) 



(44) 



where K = ~2hA. m ki J k j m£l and X = —AXE r x m {x — 
x m )R+ rn /h£lL 2 . The Pauli matrices <x operate in the 
subspace { r Di, r D 2 }. Noting that \X, K] = iC with 

C = 16VXEr /2 R%lA m /h 2 n 2 k L L, one can define the op- 



erators a = 



nc 



that satisfy [a, a t ] 



The eigen-cquations of Hi 



rzc 

, then become 



'2C 



-(at 
i(ot - 



a) 



,t 



-a) 
a 



where u n = T>\ ■ W n ,v n 
gives the equations 

at \ / u n H 
a / l u n - 



v 2 -w n 



— E n 
Combining u 



- v2C 



iv r , 



= E„ 



iVr, 



1. 



(45) 



.(46) 



Squaring Eq. (46) yields harmonic oscillator Hamilto- 
nian, and allows one to read off 



E n 



±V2Cy™ (n = 0,l,2,...). 



(47) 



Not only is there a zero-energy mode E — (the Majo- 
rana mode), but there is a ladder of localized fcrmionic 
modes, whose energy spacing is proportional to A 1 / 4 , and 
whose wavefunction components are excited harmonic os- 
cillator states. For a homogeneous gas where A = 0, the 
energy spacing becomes zero. 

In Fig. 8, we plot E n j sfnE r as a function of A 1 / 4 
(black thick curve) based on Eq. (47), and compare to 
the numerical results calculated from Eq. (22). The dot- 
dashed (red), dashed (green), and dotted (blue) curves 
show the energy levels of the first three excited states. 
We see the analytic results agree well with the numerics 
for small A. For larger A, the corrections to Eq. (44) are 
important, and the discrepancy between the analytic and 
numerical results becomes notable, especially for larger 
n. 

At n = {Eq = 0), the zero-energy mode has wave- 
function 



u{x) 
v{x) 



2ct^F 
1 

2i<jypK 



(48) 
(49) 
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FIG. 8: (color online) The gap between the MF state and 
excited states as a function of trap stiffness A 1//4 : the trap- 
ping potential is V(x) = \{x/L) 2 E r . The black (thick) curve 
is plotted based on the analytic Eq. (47). The red (dot- 
dashed), green (dashed), blue (dotted) curves are the energy 
levels E 1 /E r ,E 2 /V2E r ,E 3 /V3E r respectively. They are nu- 
merically calculated from Eq. (22) with the parameters iden- 
tical to the thick curve. 



where the width a = y A rn LE r j 'i? +m /cLA 1 / 2 is propor- 
tional to A" 1 / 4 . 

IV. SUMMARY 

We have investigated a (pseudo) spin- 1/2 spin-orbit 
(SO) coupled Fermi gas in a one-dimensional geometry. 
We first relate this system to a one-band model with 
p-wave interactions. We then described the band struc- 
ture and calculated the Berry phase 7 of the full two- 
band model. We found 7 distinguishes two topologically 
distinct sectors, with 7 = tt corresponding to a conven- 
tional superconductor. By self-consistently solving the 
Bogoliubov-de Gennes equations and calculating both 
the position resolved and momentum resolved density of 
states, we visualized the Majorana fcrmion (MF) states in 
real and momentum space at finite temperatures. These 
spectra can be probed using the position resolved or mo- 
mentum resolved radio-frequency spectroscopy [12, 18]. 
We introduced MF operators and constructed the local- 
ized MF states. We further linearized the trap near the 
location of a MF, finding an analytic expression for the 
localized MF wavefunction and the gap between the MF 
state and other edge states. 

This physics can be experimentally studied in a bundle 



of weakly coupled tubes containing fcrmionic atoms [17]. 
By applying appropriate Raman lasers to these quasi-lD 
tubes [11, 12], one can produce an array of quasi-lD SO 
coupled Fermi clouds. Our calculations show that the 
MFs can be observed in such settings. 

There are, however, significant experimental chal- 
lenges. Most notably, the Raman induced SO coupling 
relies on the ability of optical photons to flip the atomic 
hypcrfine spin. As Spiclman argues [22], if the Raman 
lasers arc detuned by a frequency A from an excited state 
multiplet (and HA is large compared to the fine struc- 
ture splitting Af), then the coupling strength ft scales as 
1/A 2 . (This is contrasted with typical AC stark shifts, 
which instead scale as 1/A. The extra suppression is due 
to quantum interference between the amplitudes arising 
from different intermediate states.) The rate of inelas- 
tic light scattering Ti also scales as 1/A 2 . The ratio 
v = r,/0 is therefore roughly independent of detun- 
ing. In terms of microscopic parameters, v oc h/AfT, 
where r is the lifetime of the excited states. For 6 Li, 
h/AfT ~ 270, for 40 K, h/AfT ~ 4.6 x 10 4 and for 87 Rb, 
h/AfT ~ 1.9 x 10 5 . One sees 40 K has a much longer life- 
time than 6 Li in a SO coupled Fermi experiment. The 
situation is even less favorable at the typical magnetic 
field ~ 830G [23] where one encounters Fcshbach reso- 
nances in 6 Li. The large magnetic field decouples the 
electron spin and the nuclear spin, and the relevant hy- 
pcrfine states effectively only differ by their nuclear spin. 
As a result, the Raman laser couplings vanish between 
these states. However for 40 K, the typical resonance field 
~ 200G [24] is much smaller, and the relevant hypcrfine 
states have larger Raman couplings. We therefore expect 

K is a promising candidate for producing an interact- 
ing SO coupled Fermi gas. 
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VI. APPENDIX 

In this Appendix, we explore the convergence of our 
self-consistent calculations with the grid spacing Sx = 



2i/n g rid- We show how the energy and the order param- 
eter for a zero-temperature homogeneous gas in a box of 
size 2L with periodic boundary conditions depends on 

"grid ■ 

Within our mean-field theory, the energy of this homo- 
geneous gas is 

E a = J2 ('o(k) UE+(k)+E-(k))) - (50) 
^~ V 1 J 9id 

where g 1D = g 1D /2L,e (k) = 4^- - //, and 
E±(k) is the excitation spectrum given in Eq. (11). 
The summation index k is discretized as fc = 

-^-(^K--,(^M^h and 
Eq. (50) can be calculated numerically. 



E g /E r 




FIG. 9: (color online) Ground state energy E g /E r versus or- 
der parameter A/E r . The blue (dashed), black (thick), red 
(dotted), green (dot-dashed) curves correspond to n gr id = 
400, 600, 800, 1000 respectively. Other parameters are gm = 
-0.02£ r , W. = 2E r , X = 0, k L L = 100, ft = E r . 

Fig. 9 shows Eg as a function of A for n gr id = 
400,600,800,1000. We find non-trivial behavior at in- 
termediate Jigrid- In particular, for these parameters and 
"■grid = 600, the energy has two local minima, and the gap 
equations has four solutions, corresponding A = and 
other three stationary points. Such behavior is an arti- 
fact of the discretization, as it goes away for n gr id > 800. 
It does, however, indicate that in the presence of an ap- 
propriate tuned optical lattice, there will be metastable 
supcrfluid states. 

In Fig. 10, we show how the order parameter A de- 
pends on n gr id- We calculate A by minimizing E g , 



We see A converges to a finite value as ngrid — > oo. For 
the simulation size n gr id = 1200 used in the main text, 
the finite grid error is |A(n KTid =oo)-A(n Kr , d =i200)| < m 

° A(n grid =oo) — 




FIG. 10: (color online) Order parameter A/E r versus 
10 3 /n gr id. The red dots are calculated from Eq. (51). The 
blue (thick) curve is an extrapolation. The parameters here 
are identical to those in Fig. 4(a) except for A = 0. 



